/*
 * ReactionSystem.cpp
 *
 *  Created on: 23 Jul 2011
 *      Author: Allan
 */

#include "ReactionSystem.h"

// C++ includes
#include <cmath>
#include <iomanip>
#include <set>

// Boost includes
#include <boost/foreach.hpp>

// Eigen includes
#include <Eigen/Core>
#include <Eigen/LU>
using namespace Eigen;

// GeoReact includes
#include "Utils/Algorithms.h"
#include "Utils/Assert.h"

ReactionSystem::ReactionSystem()
{}

ReactionSystem::ReactionSystem(const vector<Reaction>& reactions) :
reactions(reactions)
{
	// Collect all species and their indexes from the reactions
	set<tuple<string,Index>> pairs;
	
	BOOST_FOREACH(const Reaction& reaction, reactions)
		for(unsigned i = 0; i < reaction.GetSpecies().size(); ++i)
			pairs.insert(make_tuple(reaction.GetSpecies()[i], reaction.GetIndexes()[i]));
	
	species.reserve(pairs.size());
	ispecies.reserve(pairs.size());
	
	BOOST_FOREACH(auto pair, pairs)
	{
		species.push_back(get<0>(pair));
		ispecies.push_back(get<1>(pair));
	}

	// Assemble the stoichiometric matrix of the reactions
	uri = AssembleStoichiometricMatrix(species);
}

ReactionSystem::~ReactionSystem()
{}

const unsigned ReactionSystem::GetNumReactions() const
{
	return reactions.size();
}

const unsigned ReactionSystem::GetNumSpecies() const
{
	return species.size();
}

const vector<Reaction>& ReactionSystem::GetReactions() const
{
	return reactions;
}

const vector<string>& ReactionSystem::GetSpecies() const
{
	return species;
}

const vector<Index>& ReactionSystem::GetIndexes() const
{
	return ispecies;
}

const MatrixXd& ReactionSystem::GetStoichiometricMatrix() const
{
	return uri;
}

const MatrixXd ReactionSystem::AssembleStoichiometricMatrix(const vector<string>& species) const
{
	// The number of reactions (Nr) and species (Ni)
	const unsigned Nr = reactions.size();
	const unsigned Ni = species.size();
	
	MatrixXd uri = MatrixXd::Zero(Nr, Ni);
	
	for(unsigned r = 0; r < Nr; ++r) for(unsigned i = 0; i < Ni; ++i)
		uri(r, i) = reactions[r].Stoichiometry(species[i]);
	
	return uri;
}

const VectorXd ReactionSystem::EquilibriumConstants(double T, double P) const
{
	VectorXd Kr(reactions.size());
	
	for(int r = 0; r < Kr.rows(); ++r)
		Kr[r] = reactions[r].EquilibriumConst(T, P);
	
	return Kr;
}

const VectorXd ReactionSystem::ReactionQuotients(const VectorXd& a) const
{
	const unsigned Ni = species.size();
	
	// Calculate the log of the activities of the species in the system of reactions
	VectorXd ln_ai(Ni);
	
	for(unsigned i = 0; i < Ni; ++i)
		ln_ai[i] = log(a[ispecies[i]]);
	
	// Calculate the reaction quotients as Qr = exp([uri] * ln_ai)
	return (uri * ln_ai).array().exp();
}

const VectorXd ReactionSystem::KineticRates(double T, double P, const VectorXd& n, const VectorXd& a) const
{
	VectorXd Rr(reactions.size());
	
	for(int r = 0; r < Rr.rows(); ++r)
		Rr[r] = reactions[r].KineticRate(T, P, n, a);
	
	return Rr;
}

const Reaction& ReactionSystem::operator[](unsigned idx) const
{
	return reactions[idx];
}

ostream& operator<<(ostream& out, const ReactionSystem& reactions)
{
	BOOST_FOREACH(auto reaction, reactions.GetReactions())
		out << reaction << endl;
	
	return out;
}
